Liebe Studierende,
an vier Terminen in diesem Semester möchten wir Ihnen die wichtigsten Grundlagen zum Arbeiten mit dem Programmier- und Statistikprogramm R vermitteln.
Um Ihnen einen Leitfaden durch die Themen und das Semester zu geben, haben wir dieses online-Buch erstellt. Jeder Lerneinheit ist ein Kapitel zugeordnet, in dem Sie verschiedene Materialien wie Texte, Abbildungen und Videos finden.
Falls Sie Fragen haben, schreiben Sie mir gerne!
Katja Schiffers (katja.schiffers@uni-bonn.de)
aus der AG Gartenbauwissenschaften der Uni Bonn.
| Datum | Thema |
|---|---|
| 10. November | Einführung R |
| 17. November | Daten einlesen, Prinzip Anova |
| 22. Dezember | Zähldaten und Boniturdaten |
| 12. Januar | Zeitreihen und Posthoc Tests |
| 2. Februar | Bayesische Analysen |
In den nächsten Kapiteln werde ich mit Ihnen die Grundzüge der Datenauswertung mit der Statistik-Software R besprechen. Am Ende dieses Kapitels wissen Sie (hoffentlich)
R ist eine Programmiersprache, die zur Datenanalyse entwickelt wurde und damit eine Alternative zu Programmen wie SPSS, SAS und den Statistik-Funktionen von Excel. Die zwei größten Vorteile von R gegenüber point-and-click Programmen sind:
Weitere Vorteile sind:
R ist auch eine Investition in Ihre (berufliche) Zukunft. Außer Daten auswerten kann man damit
Ganz allgemein gilt R zunehmend als die Standardsprache für statistische Problemstellungen sowohl in der Wirtschaft als auch in der Wissenschaft (Amirtha et al., 2014).
Sie gelangen hier auf eine Seite mit einem Link zum R-download. Folgen Sie diesem Link und suchen Sie dann einen CRAN-mirror aus. Im Prinzip ist egal, welchen Sie wählen (auf allen sind die gleichen Ressources verfügbar, deshalb ‘Spiegel’), es macht aber Sinn, einen Mirror in der Nähe zu wählen. Es wird dann ein Installationsprogramm (z.B. eine exe-Datei oder eine .deb Datei) heruntergeladen. Wenn Sie diese öffnen, werden Sie durch die Installation geführt - bei den Auswahlmöglichkeiten passen normalerweise voreingestellten Optionen.
https://www.rstudio.com/products/rstudio/download/#download
Wenn Sie R-Studio öffnen, sehen Sie vier Fenster: den eigentlich Editor, in dem Sie R-code schreiben, die R Konsole, in der die aktive R-Session läuft, eine Übersicht über den Arbeitsspeicher und unten rechts ein Fenster in dem je nach Aktion Plots, Hilfedateien oder die Ordnerstruktur und Dateien des Verzeichnisses, in dem der R-Prozess gestartet wurde, zu sehen sind. Eventuell fehlt auch der Editor oben links noch, das ändert sich aber, wenn Sie ein neues Skript anlegen (siehe ‘Erste Schritte mit R’)
5 + 7
schreiben → auf den ‘Run’ Button klicken (über dem Editor) → der Code wird an den laufenden R Prozess in der Konsole geschickt und das Ergebnis in der Konsole ausgegeben
Ergebnis1 <- 5 + 7
eingeben → Run-Button anklicken → nichts passiert außer dass die Zeile in der Konsole auftaucht? → das Ergebnis wird in der Variable Ergebnis1 gespeichert, diese Variable wird jetzt auch in der Übersicht des Arbeitsspeichers gelistet. Wenn man in der Konsole
Ergebnis1
eingibt, wird der Wert der Variable ausgegeben. Anstelle des Pfeils kann auch ein = verwendet werden.
# Hier benutzen wir R als Taschenrechner
Ergebnis1 <- 5 + 7
kirschen <- read.csv("Kirschen.csv") eingeben um die
Daten zu importieren.Damit Sie erstmal keine Daten eingeben müssen, können Sie die
Kirsch-Daten hier herunterladen und dann mit Punkt 3 weitermachen:
https://ecampus.uni-bonn.de/goto_ecampus_file_3525090_download.html
Wenn Sie die Daten in R eingelesen haben, sollten Sie als erstes
kontrollieren, ob der Import geklappt hat. Eine erste Übersicht bekommen
Sie mit den Funktionen summary() und head().
Der Datensatz heißt ‘kirschen’, weil wir das beim Import so angegeben
haben.
summary(kirschen)
# gibt eine Zusammenfassung der Daten mit Mittelwerten etc. aus
head(kirschen)
# Zeigt die ersten 6 Zeilen des Datensatzes, wie er in R eingelesen wurde
Hier prüft man unter anderem, ob die Daten den richtigen Typ haben, also Zahlen mit denen wir rechnen wollen ‘numeric’ sind (dann werden ein paar Statistiken, wie Minimum, Mittelwert etc angegeben) und Kürzel für Versuchsbehandlungen Faktoren (dann sieht man, wie viele Beobachtungen es pro Gruppe = Level gibt). Ist das nicht der Fall, kann es jetzt geändert werden. Hier sind die Untersuchungsgruppen zum Beispiel als ‘character’ eingelesen worden (reine Zeichenketten)deshalb wandeln wir sie mit der Funktion ‘as.factor’ in Faktoren um:
kirschen$Unterlagenklon <- as.factor(kirschen$Unterlagenklon)
kirschen$Inokulat <- as.factor(kirschen$Inokulat)
Auf die einzelnen Variablen können Sie zugreifen, indem Sie zuerst den Namen des Datensatzes (= ‘kirschen’) verwenden und dann nach dem Dollar Zeichen die Spaltenüberschrift der Variable.
Auch einfache Plots sind ein gutes Mittel, um eine erste Vorstellung der Daten zu bekommen und Auffälligkeiten wie extreme Werte (verrutscher Punkt in den Zahlen) zu entdecken. Dafür eignen sich
die einfach plot-Funktion
plot(kirschen$Inokulat, kirschen$Trockengewicht,
main = "Inokulat", ylab = "Trockengewicht der Blätter [g]")oder eine schönere Version mit dem Paket ‘ggplot2’. Dazu muss das Paket ‘ggplot2’ erst installiert werden. In R-Studio geht das über ‘Tools -> Install Packages’ -> a search in the ‘Packages’ field.
library(ggplot2)
ggplot(kirschen, aes(x=Inokulat, y=Trockengewicht)) +
geom_boxplot(notch=TRUE)Hilfe zu finden ist zu Beginn des Arbeitens mit R die wichtigste Kompetenz überhaupt. Am einfachsten ist es natürlich, wenn Sie direkte Unterstützung durch jemanden bekommen, der sich mit R auskennt (und das ist bei der Bachelor-Arbeit - zumindest im Gartenbau - natürlich der Fall). Es gibt aber auch eine Reihe anderer Möglichkeiten:
?sd() → es öffnet
sich eine Hilfedatei,help.start() → öffnet Hilfeseite mit wichtigen
DokumentenNachdem Sie nun R und R-Studio kennengelernt und ein paar erste Versuche gemacht haben, stellen wir Ihnen in diesem Kapitel die wichtigsten Arbeitsschritte für ein Datenprojekt mit R vor, die Ihnen helfen, ihre Arbeit zu organisieren und die Daten in ein auswertbares Format zu bringen. Am Ende des Kapitels wissen Sie
R-Projekte sind gewissermaßen das zu Hause von allen R-Skripten und Daten, die zu einem Projekt (zum Beispiel einem Experiment) gehören. Dabei ist ein Projekt aber mehr als nur ein Ordner. Die Vorteile eines Projekts im Vergleich zu einer einzelnen neuen R Script-Datei sind, dass
kirschen <- read.csv("Kirschen.csv") immer
funktionieren, wenn die Datei ‘Kirschen.csv’ im gleicher Ordner, wie die
.Rproj-Datei liegt.Sie können ein RStudio-Projekt sowohl in einem neuen Verzeichnis als auch in einem vorhandenen Verzeichnis, in dem Sie bereits R-Code und Daten haben, erstellen. Um ein neues Projekt zu erstellen, verwenden Sie den Befehl ‘File -> New Project’.
Es öffnet sich dann ein Fenster, in dem Sie entscheiden können, ob Sie
das Projekt in einem neuen Verzeichnis, einem schon existierenden
Verzeichnis oder mit Versionskontrolle (zum Beispiel mit git) anlegen
wollen. Letzteres ist sehr sinnvoll, würde hier aber den Rahmen
sprengen.
Durch das Anlegen eines neuen Projekts in RStudio wird eine Projektdatei (mit der Erweiterung .Rproj) im Projektverzeichnis erstellt. Diese Datei kann später als Verknüpfung zum direkten Öffnen des Projekts verwendet werden.
Das Einlesen der Daten sollte jetzt kein Problem mehr sein, allerdings haben die Daten manchmal nicht den Datentyp der für bestimmte statistische Auswertungen oder Abbildungen gebraucht wird. Zuerst eine Übersicht über die gängigsten Typen:
| Bezeichnung | Beispiele |
|---|---|
| integer | 1, -2, 3, (NA) |
| numeric | 32.5, 74.3, (NA) |
| character | ‘ein Wort’, “ein Wort”, (NA) |
| logical | TRUE, FALSE, T, F, (NA) |
| factor | treat1, treat2, treat3, (NA) |
| date | 2022-10-24, (NA) |
Mit der Funktion str()kann man testen, als welche
Datentypen R die Daten interpretiert und bekommt eine Übersicht der
Struktur des Datensatzes:
str(kirschen)
Wie oben erwähnt, werden für einige Auswertungen/Abbildungen manchmal
andere Datentypen benötigt, als R erkannt hat. So werden zum Beispiel
manchmal die Bezeichnungen für Behandlungsgruppen als ‘character’
eingelesen, für eine Varianzanalyse werden aber Faktoren benötigt. In
diesem Fall kann der Datentyp mit der Funktion as.factor()
umgewandelt werden.
kirschen$Inokulat <- as.factor(kirschen$Inokulat)
Das geht natürlich nur, wenn der neue Datentyp zu der umzuwandelnden
Variablen passt: ‘Gi_12’ kann nicht in eine numerische Variable
umgewandelt werden, wohl aber die Zahlen 1, 2, 3, .. in Faktoren (die
dann als “1”, “2”, “3” ausgegeben werden). Andere häufig genutzte
Funktionen sind as.numeric() , as.character()
und as.Date().
Auch NAs können manchmal zu Fehlermeldungen oder unvollständigen Abbildungen führen. Da NA nicht Null bedeutet, sondern eher ‘wissen wir nicht’, kann zum Beispiel nicht ohne Weiteres ein Mittelwert berechnet werden.
Daten <- c(5, 3, 9, NA, 4, NA)
mean(Daten)
Die meisten Funktionen haben aber die Option, NA-Werte bei der Berechnung auszuschließen:
mean(Daten, na.rm = TRUE)
Schauen Sie also am besten zuerst in der Hilfe zu der Funktion , die
Sie verwenden möchten (?mean), ob es ein Funktionsargument
gibt, das den Umgang mit NA-Werten festlegt und ändern Sie die
Einstellung entsprechend. Wenn es gar nicht anders geht, kann man als
Alternative auch die Zeilen in den NAs vorkommen vor der Analyse
herausfiltern. Wie solche Daten-Operationen funktionieren, lernen Sie im
folgenden Kapitel.
Manchmal ist es nötig, den eingelesenen Datensatz umzustrukturieren, Teildatensätze herauszufiltern oder nue Variablen zu berechnen. Dazu sollte man am besten das Paketbündel ‘tidyverse’ installieren. Hat man das getan, steht auch eine weitere Konvertierungsfunktion zur Verfügung, die aus einem normalen Datensatz (data.frame) ein ‘tibble’ macht. Tibble haben ein paar Vorteile, zum Beispiel, dass standardmäßig nur die ersten 10 Zeilen ausgegeben werden und zu jeder Spalte der Datentyp gleich dabei steht (wobei Fließkommazahlen hier nicht als numeric sondern double (dbl) bezeichnet werden. ‘fct’ steht für Faktor, ‘int’ für Integer).
kirschen <- as_tibble(kirschen)
kirschen
Einer der wichtigsten Operatoren des tidyverse ist der Pipe-Operator
%>%. Mit ihm leiten wir ein Objekt (zum Beispiel einen
Datensatz) an eine Funktion weiter. Lesen können wir ihn als ‘dann’.
Anstelle von
Daten <- c(5, 9, 2, 5, 3, 9, 5, 3, 5)
Ergebnis <- round(sqrt(sum(Daten^2)), 3)
Ergebnis
oder
Daten_quad <- Daten^2
Daten_quad_sum <- sum(Daten_quad)
Daten_quad_sum_wurzel <- sqrt(Daten_quad_sum)
Ergebnis <- round(Daten_quad_sum_wurzel, 3)
können wir jetzt schreiben
Ergebnis <- Daten^2 %>%
sum() %>%
sqrt() %>%
round(3)
Gelesen: wir nehmen die quadrierten Daten, dann berechnen wir die Summe, dann nehmen wir die Quadratwurzel und dann runden wir auf 3 Nachkommastellen. Wir wollen hier aber keine mathematischen Berechnungen durchführen sondern Datensätze bearbeiten und umordnen. Die dafür bereitgestellten Funktionen können aber genauso mit dem Pipe-Operator genutzt werden. Hier wiederum nur die wichtisten:
| Funktion | Verwendung |
|---|---|
| drop_na() | löscht alle Zeilen des Datensatzes die NAs enthalten |
| select() | wählt Spalten aus |
| filter() | wählt Zeilen aus |
| mutate() | erstellt neue Variablen oder ersetzt scho vorhandene |
| group_by () | Berechnungen an Untergruppen von Daten |
| arrange() | sortiert den Datensatz nach einer bestimmten Variable |
| pivot_wider() | verringert die Anzahl der Zeilen, erhöht die Anzahl der Spalten |
| pivot_longer() | erhöht die Anzahl der Zeilen, verringert die Anzahl der Spalten |
Der Funktionsname spricht für sich: alle Zeilen des Datensatzes, die ein NA enthalten werden komplett aus dem Datensatz gelöscht:
library(tidyverse)
# mit der Funktion 'tibble()' legt man einen Datensatz an
Daten1 <- tibble(spalte1 = c(1, NA, 5, 7), spalte2 = c("a", "b", NA, "d"))
# wir schauen uns den Datensatz an:
Daten1
# wir löschen die Zeilen mit NA heraus
Daten2 <- Daten1 %>% drop_na()
# uns sehen uns das Ergebnis an:
Daten2
Mit der Funktion select() können sehr einfach Spalten ausgewählt werden. Hier verwenden wir zusätzlich die Funktion head() um nur die ersten 6 Spalten ausgeben zu lassen:
kirschen %>% select(Trockengewicht, Blattfläche) %>% head()
Mit einem Minus vor der Spaltenüberschrift, werden die entsprechenden Spalten ausgeschlossen:
kirschen %>% select(-Trockengewicht, -Blattfläche) %>% head()
Wenn der Teildatensatz gespeichert werden soll, muss er einem Variablennamen zugewiesen werden
kirschen_neu <- kirschen %>% select(-Trockengewicht, -Blattfläche)
head(kirschen_neu)
…funktioniert analog mit Zeilen. Hier wird allerdings mit einem doppelten Gleichzeichen geprüft, welche Bedingung erfüllt sein muss, um im Datensatz zu bleiben.
kirschen_nurPDV <- kirschen %>% filter(Inokulat == "PDV")
summary(kirschen_nurPDV)
Aber Achtung! Wie man in der summary sieht, sind trotzdem noch alle 4 Level von ‘Inokulat’ vorhanden, nur sind die Zähler der Level ‘ohne’, ‘PDV + PNRSV’ und ‘PNRSV’ auf 0 gesetzt. In Analysen und bei Abbildungen würden sie deshalb trotzdem auftauchen. Deshalb ist es wichtig, die nicht mehr benötigten Level mit der Funktion droplevels() herauszuwerfen:
plot(kirschen_nurPDV$Inokulat, kirschen_nurPDV$Frischgewicht)
kirschen_nurPDV <- kirschen %>%
filter(Inokulat == "PDV") %>%
droplevels()
plot(kirschen_nurPDV$Inokulat, kirschen_nurPDV$Frischgewicht)
Genau entgegengesetzt zum ==funktioniert das ‘ungleich’
Symbol !=. Hiermit würden alle Zeilen im Datensatz bleiben
die NICHT “PDV” als Inokulat-Level haben.
Mit mutate() kann man eine neue Variable anlegen oder eine schon vorhandene überschreiben. Das Beispiel erklärt am besten, wie es funktioniert:
kirschen %>% mutate(Gesamttrieblänge_m = Gesamttrieblänge/100)
… wird dazu verwendet, aus einem Datensatz im ‘weiten’ Format ein Datensatz im ‘langen’ Format zu machen. Häufig ist es zum Beispiel so, dass es bei Zeitreihen (Wachstum eines Baumes) für jeden Messzeitpunkt eine eigene Spalte gibt in der die Größe eingetragen wird. Manchmal ist es aber notwendig, dass alle gemessenen Größen in einer einzigen Spalte untereinander stehen und in der Spalte daneben, der Zeitpunkt als Variable angegeben ist. Wird im Beispiel klarer…
Wachstum <- tibble(Baum = c("Nr1", "Nr2", "Nr3", "Nr4"),
Zeitp1 = c(338,459,558,439),
Zeitp2 = c(437,560,729,658),
Zeitp3 = c(469,579,937,774))
Wachstum_lang <- Wachstum %>%
pivot_longer(
cols = c(Zeitp1, Zeitp2, Zeitp3),
names_to = "Zeitpunkt",
values_to = "Größe_cm")
# cols: die Zeilen die zusammengefügt werden sollen
# names_to: neuer Spaltenname des Faktors
# values_to: neuer Spaltenname der Daten
# 'Zeitpunkt' und 'Baum' wurden noch nicht als Faktor erkannt, deshalb:
Wachstum_lang$Zeitpunkt <- as.factor(Wachstum_lang$Zeitpunkt)
Wachstum_lang$Baum <- as.factor(Wachstum_lang$Baum)
Dieses lange Format wird zum Beispiel benötigt, um einen ggplot, der das Wachstum darstellt, machen zu können
ggplot(Wachstum_lang, aes(x = Zeitpunkt, y = Größe_cm, color=Baum)) +
geom_point() + labs(y="Größe (cm)", x="Zeitpunkt")
Und ganz zum Schluss noch die Umkehrfunktion von pivot_longer(): pivot_wider. Sie nimmt einen Faktor und eine Messvariable und “verteilt” die Werte der Messvariable über neue Spalten, welche die Stufen des Faktors repräsentieren:
Wachstum_ww <- Wachstum_lang %>%
pivot_wider(
names_from = "Zeitpunkt",
values_from = "Größe_cm"
)
Bevor wir in die eigentliche Datenanalyse einsteigen, ist es wichtig, dass Sie das Grundprinzip der statistischen Tests verstanden haben. Auf die Gefahr hin, dass es eine Wiederholung für einige von Ihnen ist, wollen wir deshalb nochmal die Idee der Varianzanalyse, die den meisten Tests zugrunde liegt, erklären. Am Ende dieses Kapitels
Wir nehmen an, wir haben einen sehr einfachen Versuch mit Tomaten durchgeführt, in dem die Hälfte aller 40 Pflanzen gedüngt wird und die andere Hälfte als Kontrollgruppe ungedüngt bleibt. Am Ende bestimmen wir das Gesamtgewicht der Tomaten pro Pflanze als ein Maß für den Ertrag. Jetzt interessiert uns, ob die Düngung einen signifikanten Effekt auf den Ertrag hat. Unsere Hypothese ist: Die Düngung hat einen Effekt auf den Ertrag. Statistisch ausgedrückt bedeutet das, dass die Ertragsdaten aus zwei Normalverteilungen mit unterschiedlichen Mittelwerten stammen. Diese Hypothese ist auf folgender Abbildung auf der linken Seite dargestellt. Das zugehörige statistische Modell hat zwei Parameter: a ist der Mittelwert der Normalverteilung ohne Dünger, b die Differenz zwischen a und dem Mittelwert für die Normalverteilung mit Dünger (die Effektgröße). Die Variable Dünger nimmt entweder den Wert 0 (ohne Dünger) oder 1 (mit Dünger) an. Die Streuung der Daten (sie liegen nicht alle genau auf dem Mittelwert) wird durch den Fehlerterm aufgefangen.
Die Null-Hypothese zu diesem Modell ist: Die Düngung hat keinen Effekt auf den Ertrag -> die Daten stammen alle aus einer einzigen Normalverteilung mit dem Mittelwert a (rechte Seite der Abbildung). Entsprechend einfacher ist das statistische Modell.
In R formulieren wir diese beiden konkurrierenden Modelle mit der
Funktion lm(). lm steht für ‘linear model’ und es schätzt
Mittelwert (oder in anderen Fällen auch Achsenabschnitt und Steigung
einer Regressionsgeraden), indem es die Werte so wählt, dass die
Fehlerquadrate (quadrierter Abstand der Messungen zu den Mittelwerten)
möglichst klein sind.
# hier lesen wir keine Daten aus Excel ein sondern geben sie einfach händisch direkt
# in R ein. Zuerst Ertragsdaten von 10 ungedüngten Pflanzen
ohne <- c(417.1,558.8,519.1,611.7,450.3,461.9,517.3,453.2,533.1,514.7)
# dann von 10 gedüngten Pflanzen
mit <- c(581.9,517.4,641.3,659.8,587.8,483.3,703.0,589.5,532.4,569.3)
# hier werden die beiden Datensätze mit c() wie combine
# zu einem Vektor (= Spalte) zusammengeführt
Ertrag <- c(ohne, mit)
# jetzt generieren wir noch einen Vektor der die Behandlung anzeigt
Düngung <- gl(2, 10, 20, labels = c("ohne","mit"))
# Mit diesen beiden Vektoren können wir das statistische Modell formulieren,
# dass der Ertrag von der Düngung abhängt. In R wird das mit der Tilde ~ ausgedrückt
# (oben links auf der Tastatur).
# Wir nennen dieses Modell 'model1'
model1 <- lm(Ertrag ~ Düngung)
# mit der folgenden Zeile schauen wir uns eine summary (Zusammenfassung)
# des gefitteten Modells ansummary(model1)
In der Zeile ‘Intercept’ finden wir unter ‘Estimate’ den geschätzten Mittelwert für die Gruppe ohne Düngung, also das a in der Abbildung oben, linke Seite. Darunter, in der Zeile ‘Düngungmit’ ist die Differenz (das b aus der Abbildung) das auf das a aufaddiert werden muss, um die Schätzung für den Mittelwert der gedüngten Pflanzen zu erhalten.
Als nächstes formulieren wir die Nullhypothese
##
## Call:
## lm(formula = Ertrag ~ 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -128.04 -38.30 -12.39 43.08 157.85
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 545.15 16.68 32.69 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 74.59 on 19 degrees of freedom
Wie erwartet, erhalten wir jetzt nur einen Mittelwert über alle Tomatenpflanzen.
Welches der beiden Modelle ist nun besser? Das Modell mit Düngung als erklärender Variable passt sich besser an die Daten an -> kleinere Summe der Fehlerquadrate, ist aber komplizierter.
Um zu entscheiden, ob eine erklärende Variable einen signifikanten
Effekt hat, stellen wir die Frage: Wie wahrscheinlich ist es, die Daten
so - oder mit noch größeren Unterschieden zwischen den Behandlungen - zu
beobachten, wenn in Wahrheit die erklärende Variable keinen Effekt hat
(= die Nullhypothese wahr ist bzw. alle Daten aus einer einzigen
Normalverteilung stammen). Um diese Wahrscheinlichkeit ausrechnen zu
können, müssen die Daten innerhalb der Behandlungsgruppen normalverteilt
sein und die Streuung (Varianz) muss in etwa gleich sein. Das sind die
wichtigsten Vorraussetzungen für eine Varianzanalyse, weil sonst evtl.
die berechnete Wahrscheinlichkeit nicht stimmt. In R nutzten wir die
Funktion anova() um die beiden Modelle zu vergleichen und
die oben genannte Wahrscheinlichkeit auszurechnen:
anova(model1, null_model)
Die Wahrscheinlichkeit nach der wir suchen, finden wir unter Pr(>F) und ist der berühmte p-Wert. Die Wahrscheinlichkeit zufällig solche (oder noch extremere) Daten zu finden, wie wir sie für die Tomaten gefunden haben, obwohl die Düngung in Wahrheit keinen Einfluss auf den Ertrag hat, ist nur 0,86% und damit so unwahrscheinlich, dass wir den Effekt der Düngung als eindeutig signifikant interpretieren können. Damit ist model1 das bessere Model. Eine allgemeine, wenn auch etwas arbiträre Übereinkunft, ist, dass alle Variablen mit einem p-Wert von unter 5% bzw p < 0.05 als signifikant gelten.
Kritik an der Verwendung des p-Wertes
In den letzten Jahren gibt es mehr und mehr Kritik an dem starken Fokus auf den p-Wert. Zum einen, weil zumeist einfach das Signifikanzniveau verwendet wird, was fast immer verwendet wird, ohne weiter darüber nachzudenken: 5%. Dieses Niveau stellt einen Kompromiss zwischen der Fähigkeit eine Entdeckung zu machen und unserer Bereitschaft ein fehlerhaftes Testergebnis zu akzeptieren dar. Es sollte eigentlich einleuchtend sein, dass das Signifikanzniveau für verschiedene Fragestellungen auch unterschiedlich angesetzt werden sollte (was der Begründer der frequentistischen Statistik, Ronald Fisher auch schon 1956 selbst so geraten hat). Zum anderen ist der Nachweis eines signifikanten Effekts außer von der Effektgröße stark vom Informationsgehalt der Daten abhängig (wie groß ist die Stichproble). So kann es sein, dass, wenn man hunderte von Proben hat, eine Variable als signifikant nachgewiesen werden kann, ihr Effekt und die wissenschaftliche Relevanz aber nur minimal ist. Es ist also immer, immer wichtig, den Bezug zur eigentlichen Frage und zum untersuchten System zu behalten und die Ergebnisse mit gesundem Menschen- und Sachverstand zu interpretieren.
Wie oben schon kurz angesprochen, müssen für glaubwürdige Ergebnisse einer Varianzanalyse einige Vorraussetzungen erfüllt sein
Um die Normalverteilung der Daten innerhalb der Behandlungsgruppen zu
prüfen, ist es einfacher, sich die Daten nicht gruppenweise anzuschauen,
sondern einfach die Residuen/Fehler aller Daten (also den Abstand jedes
Datenpunktes zum geschätzten Gruppen-Mittelwert) zu nehmen. Es gibt
numerische Tests, um die Normalverteilung der Residuen zu prüfen, in den
letzten Jahren setzt sich aber mehr und mehr eine visuelle Überprüfung
durch. Am besten eignet sich dafür ein qq-Plot (Quantile-Quantile-Plot).
Wenn die Punkte annähernd auf einer Geraden liegen, sind die Residuen
normalverteilt. Da ‘annähernd’ ein dehnbarer Begriff ist, liefert die
Funktion qqPlot aus dem Package ‘car’ ein
Konfidenzintervall. Liegen die Punkte innerhalb dieses Intervalls, kann
man von einer Normalverteilung ausgehen. Das Paket ‘car’ müssen Sie in
R-Studio vorher installieren (Tools -> Install Packages) Im
jupiter-Notebook sollten Sie zu Beginn der Zelle diesen Code
einfügen
# Schnelles Repository: Posit Public Package Manager
options(repos = c(PPM = "https://packagemanager.posit.co/cran/__linux__/jammy/latest"))
# Hier die benötigten Pakete auflisten
pkgs <- c("car")
# Nur fehlende Pakete installieren
to_install <- setdiff(pkgs, rownames(installed.packages()))
if (length(to_install)) install.packages(to_install)Dann geht es weiter:
## Loading required package: carData
##
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
##
## recode
## The following object is masked from 'package:purrr':
##
## some
## [1] 17 4
Außer dem Plot gibt die Funktion die ‘Namen’ der Punkte aus, die im Bezug auf die x-Achse outlier darstellen (also Punkte, die besonders stark abweidchen,hier Punkt 4 und Punkt 17). Es lohnt sich, diese Punkte genauer anzuschauen (waren die Versuchsbedingungen hier aus irgendeinem Grund anders?) und evtl. auszuschließen, insbesondere, wenn sie außerdem außerhalb der Konfidenzlinien liegen.
Eine gute Erklärung wie qq-plots Funktionieren finden Sie hier: https://www.youtube.com/watch?v=okjYjClSjOg … auf Englisch und mit toller Intro :)
Um die Homoskedastizität visuell zu überprüfen, plotten Sie die
Residuen gegen die von der Funktion lm() geschätzten Mittelwerte. Wenn
die Daten der Gruppen etwa einen gleichen Bereich auf der y-Achse
umfassen, sind die Varianzen homogen. Hierfür gibt es leider derzeit
keine Konfidenzintervalle, eine Daumenregel ist aber, dass die Varianzen
als homogen gelten, wenn die größte Varianz nicht mehr als 1,5fach
größer ist als die kleinste Varianz. In unserem Beispiel reichen die
Residuen beider Gruppen von etwa -100 bis 100, die Varianzen sind also
annähernd gleich groß. Formal kann Varianzhomogenität mit dem Leven’s
Test überprüft werden. R bietet dazu die Funktion
leveneTest() aus dem package car. Wer
interessiert daran ist, kann sich mit ?leveneTest die
Hilfeseite mit Erklärungen zur Anwendung der Funktion anschauen.
Für unsere Analyse sind alle Annahmen erfüllt. Wir können also davon ausgehen, dass die Düngung in unserem hypothetischen Experiment tatsächlich einen signifikanten Effekt auf den Ertrag hat. In einer Arbeit würden wir den F-Wert (aus summary(model1)), die Anzahl der Freiheitsgrade und den p-Wert angeben, wobei p-Werte, die kleiner als 0.05 sind, üblicherweise mit p < 0.05 angegeben werden.
In der letzten Woche haben Sie gelernt, wie man Experimente mit kontinuierlichen Daten der abhängigen Variable auswertet. In den Agrarwissenschaften kommen aber auch häufig andere Datentypen vor. In diesem Kapitel beschäftigen wir uns speziell mit Zähldaten (wieviele Äpfel trägt der Baum? wieviele Läuse sitzen auf der Pflanze? wieviele Eier legt das Huhn?)
Am Ende dieses Kapitels
Im letzten Kapitel haben wir gesehen, dass eine Voraussetzung für die Durchführung einer ANOVA die Normalverteilung der Residuen (Fehler, Abstände der Messpunkte zum Mittelpunkt) ist. Bei kontinuierlichen Daten wie Größe, Gewicht etc. ist das normalerweise der Fall, wenn das Modell richtig formuliert ist (d.h. wenn alle Faktoren, die einen Einfluss haben, enthalten sind und deren Beziehung - additiv, multitplikativ, exponentiell - korrekt dargestellt ist). Sind die Daten nicht kontinuierlich, können wir erstmal nicht von einer Normalverteilung der Residuen ausgehen. Häufig sind auch die Varianzen (= Streuung der Daten) nicht homogen über die Behandlungsgruppen. Das ist zum Beispiel bei Zähldaten der Fall, zumindest, wenn es sich um eher kleine Zahlen handelt (bis etwa 7).
In so einem Fall wurde früher oft dazu geraten, die Daten zu transformieren, um sie ANOVA-fähig zu machen. Das birgt allerdings verschiedene Schwierigkeiten:
Die meisten Ratschläge, Daten zu transformieren, stammen aus der Zeit, als Varianzanalysen noch von Hand durchgeführt wurden und kompliziertere Modelle kaum zu rechnen waren. Zum Glück haben wir heute Computer und R, um auch viele nicht-normalverteilte Daten mittels einer ANOVA (Varianzanalyse) analysieren zu können. Der Trick hierbei ist, nicht die schon bekannten linearen Modelle (lm) zu verwenden, sondern generalisierte lineare Modelle (glm).
… können Poisson-verteilte Residuen (aus Zähldaten), binomial-verteilte Residuen (aus Anteilsdaten) und viele andere Residuenverteilungen beschreiben. Ich werde hier nicht auf die Charakteristika dieser Verteilungen eingehen, weil es für das Verständnis der Auswertung nicht essentiell ist. Wenn Sie Interesse haben, finden Sie aber gute Erklärungen über eine Internet-Suche.
In einem normalen linearen Modell (oberer Teil der Abbildung) liefern die erklärenden Variablen (zum Beispiel ‘Düngung’)und die Modellgleichung direkt eine Schätzung der abhängigen Variable \(\mu\) (zum Beispiel ‘Ertrag). Die Abstände zu den tatsächlich beobachteten Daten \(y\) sind die Fehler oder Residuen und es wird angenommen, dass sie normalverteilt sind. Bei einem GLM wird zusätzlich eine ’link-function’ benötigt, die das Ergebnis der Modellgleichung auf den biologisch sinnvollen Bereich abbildet. So wird zum Beispiel verhindert, dass negative Zahlen geschätzt werden (ein Baum mit -5 Äpfeln). Zusätzlich ist die Annahme der Verteilung der Residuen/Fehler flexiber, wie oben bereits erwähnt.
Exkurs: Schätzung der Modell-Parameter
Grundsätzlich findet man die Modell-Parameter (bei einer Geradengleichung \(y = a * x + b\), wären das zum Beispiel \(a\) und \(b\)), für die die Daten am wahrscheinlichsten sind, indem man das Maximum-Likelihood Prinzip anwendet, sprich die deviance - also die Abweichung der geschätzten von den beobachteten Werten - minimiert. Bei normalverteilten Fehlern ist die deviance als die Summe der Fehlerquadrate RSS (residual sum of squares) definiert: \(\sum (\mu - y)^2\). Bei nicht-normalverteilten Fehlern funktioniert der Trick mit dem Minimieren der Summe der Fehlerquadrate nicht. Die deviance ist hier etwas komplizierter definiert. Für Poisson-verteilte Daten zum Beispiel \(2\sum(y_i log(y/\mu)-(y-\mu))\).
Leichter verständlich, wird es sicher, wenn wir uns ein konkretes Beispiel anschauen - hier können Sie einen Datensatz zu einem Tomatenversuch herunterladen: https://ecampus.uni-bonn.de/goto_ecampus_file_2786370_download.html Speichern Sie ihn in ihrem R-Projekt (wenn Sie eins angelegt haben) im Ordner ‘data’.
# Zuerst lesen wir den Tomaten-Datensatz ein
tomaten_daten <- read.csv("data/Tomaten.csv")
# Eine kurze Kontrolle, ob alles funktioniert hat
summary(tomaten_daten)## Datum Zeitpunkt Haus Bedachung
## Length:865 Min. : 1.000 Min. :1.000 Length:865
## Class :character 1st Qu.: 3.000 1st Qu.:2.000 Class :character
## Mode :character Median : 5.000 Median :4.000 Mode :character
## Mean : 5.089 Mean :3.504
## 3rd Qu.: 7.000 3rd Qu.:5.000
## Max. :10.000 Max. :6.000
##
## Pflanzen.Nr Veredelung Salzstress Trieblänge
## Min. : 1.00 Length:865 Length:865 Min. : 0.0
## 1st Qu.:25.00 Class :character Class :character 1st Qu.: 80.0
## Median :49.00 Mode :character Mode :character Median :104.0
## Mean :48.55 Mean :106.4
## 3rd Qu.:72.00 3rd Qu.:130.0
## Max. :96.00 Max. :196.0
##
## Anzahl_Blätter Anzahl_Blütenstände Anzahl_grüne_Früchte Anzahl_rote_Früchte
## Min. : 0.00 Min. : 0.000 Min. : 0.00 Min. : 0.000
## 1st Qu.:15.00 1st Qu.: 3.000 1st Qu.: 2.00 1st Qu.: 0.000
## Median :19.00 Median : 4.000 Median : 7.00 Median : 0.000
## Mean :19.62 Mean : 4.467 Mean :11.14 Mean : 1.029
## 3rd Qu.:24.00 3rd Qu.: 6.000 3rd Qu.:17.00 3rd Qu.: 0.000
## Max. :36.00 Max. :36.000 Max. :58.00 Max. :15.000
##
## Anzahl_Blütenendfäule Anzahl_geerntete_Früchte Gewicht_Trieb Gewicht_Blätter
## Min. : 0.000 Min. : 0.000 Min. : 0.0 Min. : 0
## 1st Qu.: 0.000 1st Qu.: 0.000 1st Qu.:235.0 1st Qu.: 820
## Median : 1.000 Median : 0.000 Median :368.0 Median :1237
## Mean : 3.042 Mean : 0.341 Mean :328.6 Mean :1077
## 3rd Qu.: 5.000 3rd Qu.: 0.000 3rd Qu.:428.0 3rd Qu.:1418
## Max. :17.000 Max. :10.000 Max. :579.0 Max. :1783
## NA's :768 NA's :768
## Gewicht_grüne_Früchte Gewicht_rote_Früchte Gewicht_Blütenendfäule
## Min. : 0.0 Min. : 0.0 Min. : 0.00
## 1st Qu.: 365.0 1st Qu.: 0.0 1st Qu.: 0.00
## Median : 688.0 Median : 0.0 Median : 0.00
## Mean : 633.2 Mean : 18.9 Mean : 21.93
## 3rd Qu.: 964.0 3rd Qu.: 0.0 3rd Qu.: 28.00
## Max. :1222.0 Max. :312.0 Max. :221.00
## NA's :768 NA's :768 NA's :768
# Es fällt auf, dass die beiden erklärenden Variablen 'Veredelung' und 'Salzstress',
# die wir gleich verwenden wollen, as 'character' eingelesen worden sind.
# Das kann bei den Abbildungen zu Problemen führen. Deshalb wandeln wir sie
# zuerst in 'factors' um:
tomaten_daten$Veredelung <- as.factor(tomaten_daten$Veredelung)
tomaten_daten$Salzstress <-
as.factor(tomaten_daten$Salzstress)
# Für die folgende Abbildung benötigen wir das package 'lattice'
# (muss wie immer zuerst installiert und dann geladen werden).
library(lattice)
# Um die Daten ein bisschen kennenzulernen (und noch ein paar neue Arten von plots kennenzulernen),
# bilden wir sie erstmal in einem xyplot ab. Insbesondere interessieren wir uns
# jetzt für die Anzahl an Früchten, die Blütenendfäule haben in Abhängigkeit von
# Salzstress und Veredelung.
#
with(tomaten_daten, xyplot(Anzahl_Blütenendfäule ~ Salzstress, groups = Veredelung))Hier sehen wir nicht besonders viel, außer dass die Varianz sehr hoch ist. Machen wir noch einen Interaktionsplot, in dem die Mittelwerte der Gruppen (mit/ohne Salzstress, mit/ohne Veredelung) berechnet und abgebildet werden.
# hier machen wir ein Interaktionsplot
with(tomaten_daten, interaction.plot(x.factor = Veredelung,
trace.factor = Salzstress, response = Anzahl_Blütenendfäule))Jetzt erkennen wir, dass veredelte Tomaten mit Salzstress eine
deutlich geringere Zahl von Früchten mit Blütenendfäule haben, als ohne
Salzstress. Bei nicht-veredelten Pflanzen ist der Effekt genau
umgekehrt, allerdings deutlich kleiner. Es sieht also so aus als gäbe es
eine Interaktion zwischen den beiden erklärenden Variablen Salzstress
und Veredelung: der Effekt der einen erklärenden Variablen auf die
abhängige Variable hängt vom Wert der anderen erklärenden Variablen ab.
Um diesen Effekte auf Signifikanz zu testen, fitten wir jetzt
generalisierte lineare Modelle, die mit der Funktion glm()
aufgerufen werden:
glm_model_1 <- glm(Anzahl_Blütenendfäule ~ Veredelung * Salzstress,
family = poisson, data = tomaten_daten)Wie Sie sehen, funktioniert das ganz ähnlich wie mit der Funktion
lm(). Es muss lediglich der zusätzliche Parameter ‘family’
angegeben werden. Hier wählen wir ‘poisson’, weil Zähldaten
Poisson-verteilt sind.
Aber Vorsicht! Bevor wir das Modell weiter analysieren können müssen wir uns noch das Verhältnis von Residual deviance (die Fehler, die es nach der Schätzung der Mittelwerte noch gibt) und den Freiheitsgraden/degrees of freedom anschauen: Ein GLM nimmt eine bestimmte Beziehung zwischen der Residual deviance und der Modellvorhersage an. Bei vielen agrarwissenschaftlichen Daten ist die Varianz aber größer als unter Poisson-Verteilung erwartet (z.B. weil nicht gemessene Variablen die abhängige Variable beeinflussen) -> das bezeichnet man als Overdispersion. Wenn die Residual deviance (zweitletzte Zeile in der summary) mehr als 1,5 mal höher ist als die Freiheitsgrade, liegt Overdispersion vor.
##
## Call:
## glm(formula = Anzahl_Blütenendfäule ~ Veredelung * Salzstress,
## family = poisson, data = tomaten_daten)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.01664 0.04093 24.840 < 2e-16 ***
## Veredelungohne 0.08507 0.05675 1.499 0.134
## Salzstressohne 0.27847 0.05414 5.143 2.70e-07 ***
## Veredelungohne:Salzstressohne -0.37364 0.07854 -4.757 1.96e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 4293.7 on 864 degrees of freedom
## Residual deviance: 4256.1 on 861 degrees of freedom
## AIC: 5833.3
##
## Number of Fisher Scoring iterations: 6
4256/861 ist deutlich größer als 1,5, wir haben in den Daten also tatsächlich Overdispersion.
Lösungen:
glm_model_1 <- glm(Anzahl_Blütenendfäule ~ Veredelung * Salzstress,
family = quasipoisson, data = tomaten_daten)In diesem Modell wird nun zusätzlich geschätzt, wie hoch die Varianz in den einzelnen Behandlungsgruppen ist. Hiermit können wir weiterarbeiten: um die Interaktion zwischen Veredelung und Salzstress auf Signifikanz zu testen, formulieren wir ein einfacheres, genestetes Modell, in dem diese Interaktion fehlt (wenn sich ein Modell B von einem Modell A ableiten lässt, dass heißt, ein oder mehrere erklärende Variablen oder Interaktionen herausgenommen worden sind, sagt man, das Modell B ist in Modell A genestet). Dazu ersetzen wir das * (das für Interaktion steht) mit einem +. So haben beide erklärenden Variablen einen unabhängigen Einfluss auf die Anzahl der Früchte mit Blütenendfäule aber keine Interaktion mehr.
glm_model_2 <- glm(Anzahl_Blütenendfäule ~ Veredelung + Salzstress,
family = quasipoisson, data = tomaten_daten)Mit der Funktion anova() vergleichen Sie die Modelle und
testen, ob die Interkation signifikant ist. Zusätzlich müssen Sie hier
die Statistik angeben, mit der die Modelle verglichen werden sollen,
weil das bei glms mit eindeutig ist. Für quasipoisson-verteilte
Residuen, wird der F-test benötigt.
## Analysis of Deviance Table
##
## Model 1: Anzahl_Blütenendfäule ~ Veredelung * Salzstress
## Model 2: Anzahl_Blütenendfäule ~ Veredelung + Salzstress
## Resid. Df Resid. Dev Df Deviance F Pr(>F)
## 1 861 4256.1
## 2 862 4278.8 -1 -22.723 4.8598 0.02775 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Wie schon vermutet, ist die Interaktion signifikant. Wenn eine Interaktion signifikant ist, bleiben auch die beiden erklärenden Variablen, die interagieren auf jeden Fall im Modell. Die Varianzanalyse ist also jetzt abgeschlossen.
Wir könnten uns jetzt mit summary(*Modellname*) noch die
geschätzten Parameter-Werte anschauen. Hier müssen Sie aber beachten,
dass sie auf der link-Skala angegeben sind (siehe oben, die Parameter
werden so angepasst, dass das Modell sinnvolle Ergebnisse liefert) und
zurück-transformiert werden müssen, um den eigentlichen Wert zu
erhalten. Am einfachsten ist es, die Funktion predict() zu
verwenden: Als erstes Argument nennen wir den Namen des glm-Modells,
dann erstellen wir einen Datensatz mit den Variablen und Werten, für die
wir die Vorhersagen sehen möchten, also mit/ohne Salzstress und mit/ohne
Veredelung in allen 4 möglichen Kombinationen. Als type
wählen wir ‘response’ aus, weil wir die Schätzungen nicht auf der
link-Skala sondern auf der response-Skala (das sind die
zurücktransformierten, interpretierbaren Werte, siehe oben) haben
wollen.
predict(glm_model_1, data.frame(Salzstress = c("mit", "mit", "ohne", "ohne"),
Veredelung = c("mit", "ohne", "mit", "ohne")), type = "response")## 1 2 3 4
## 2.763889 3.009302 3.651376 2.736111
Der geschätzte Wert für die Anzahl an Früchten mit Blütenendfäule für Pflanzen mit Salzstress und mit Veredelung ist also 2,76, für Pflanzen mit Salzstress und ohne Veredelung 3.01 usw.
In dieser Woche geht es gleich um zwei Datentypen: Anteils- und Bonitur-Daten. Anteilsdaten können mit generalisierten linearen Modellen (GLMs) ausgewertet werden, die Sie beim letzten Mal schon kennengelernt haben - allerdings mit zwei kleinen Änderungen. Bonitur-Daten sind im Gartenbau ziemlich häufig, sind allerdings in Bezug auf die Auswertung ein eher undankbarer Datentyp: oft ist es schwierig bis unmöglich, sie durch Verteilungen wie der Normalverteilung zu charakterisieren. Deshalb wiedersetzen sie sich auch den Auswertungen, die auf solchen Verteilungen beruhen.
Am Ende dieses Kapitels wissen Sie
Im Keimungsversuch mit Radis und Bohnen auf zwei verschiedenen Substraten haben Sie Anteilsdaten aufgenommen: wie viele der 80 ausgesäten Samen sind gekeimt? Die csv-Dateien finden Sie hier: https://ecampus.uni-bonn.de/goto_ecampus_file_2786369_download.html
Wie Sie sich schon denken können, erfüllt auch dieser Datentyp nicht die Annahmen, die für eine ‘normale’ ANOVa gegeben sein müssen. Zum Glück können wir wieder auf die generalisierten linearen Modelle zurückgreifen, allerdings mit zwei kleinen Änderungen. Aber fangen wir von vorne an.
Zuerst laden wir das Paket ‘tidyverse’, das uns zum Beispiel das Umstrukturieren von Datensätzen erleichtert, lesen die Daten ein (bei mir sind sie wieder im Ordner ‘data’ innerhalb des Ordners, in dem sich das r-Skript befindet gespeichert) und kontrollieren, ob alles geklappt hat:
library(tidyverse)
Ansaaten <- read.csv("data/Ansaaten.csv", header = T)
head(Ansaaten)
str(Ansaaten)
summary(Ansaaten)Um uns die Keimungsraten genauer anzuschauen, erstellen wir einen neuen Datensatz ‘Keimung’, indem wir mit dem pipe-Operator %>% aus dem Datensatz ‘Ansaaten’ die entsprechenden Spalten auswählen:
## Art ID Substrat Ausfaelle Angelaufen
## 1 Bohne 1 Presstopf 5 79
## 2 Bohne 2 Presstopf 3 81
## 3 Bohne 3 Presstopf 6 78
## 4 Bohne 4 Presstopf 6 78
## 5 Bohne 5 Presstopf 21 63
## 6 Bohne 6 Presstopf 23 61
Um einen Säulendiagramm zu erstellen, in dem die Anzahl der angelaufenen Samen und der Ausfälle gestapelt dargestellt werden, müssen wir den Datensatz umstrukturieren, um für jede Anzahl ‘Ausfaelle’ oder ‘Angelaufen’ eine eigene Spalte zu bekommen. Dazu nutzen wir die Funktion ‘reshape’. Wichtig sind hierbei die Parameter
Keimung_w <- reshape(Keimung,
direction = "long",
varying = list(names(Keimung)[4:5]),
v.names = "Anzahl",
idvar = c("Art", "ID", "Substrat"),
timevar = "Beobachtung",
times = c("Ausfaelle", "Angelaufen"))
head(Keimung_w)## Art ID Substrat Beobachtung Anzahl
## Bohne.1.Presstopf.Ausfaelle Bohne 1 Presstopf Ausfaelle 5
## Bohne.2.Presstopf.Ausfaelle Bohne 2 Presstopf Ausfaelle 3
## Bohne.3.Presstopf.Ausfaelle Bohne 3 Presstopf Ausfaelle 6
## Bohne.4.Presstopf.Ausfaelle Bohne 4 Presstopf Ausfaelle 6
## Bohne.5.Presstopf.Ausfaelle Bohne 5 Presstopf Ausfaelle 21
## Bohne.6.Presstopf.Ausfaelle Bohne 6 Presstopf Ausfaelle 23
Jetzt können wir den Plot erstellen
ggplot(Keimung_w, aes(x=factor(ID), y=Anzahl, fill=Beobachtung)) +
facet_grid(. ~ Art * Substrat) +
geom_bar(position="stack", stat="identity") +
scale_fill_manual("legend",
values = c("Angelaufen" = "darkgreen", "Ausfaelle" = "lightgrey")) +
geom_col(position = position_stack(reverse = TRUE))Wir sehen, dass die Erfolgsrate bei den Radis unabhängig vom Substrat sehr hoch ist. Bei den Bohnen gibt es mehr Ausfälle, insbesondere bei den auf Steinwolle ausgesäten.
Sind diese Unterschiede signifikant? Um das zu testen, können wir wieder die generalisierten linearen Modelle fitten und mittels der anova-Funktion vergleichen. Bei Anteilsdaten wie diesen erfolgt davor allerdings noch ein Vorbereitungsschritt: Wir generieren aus den beiden Spalten ‘Angelaufen’ und ‘Ausfaelle’ eine einzige abhängige Variable, die wir ‘Keimungsrate’ nennen. Es ist wichtig, die Rate nicht selber auszurechnen: ‘1 von 10 Samen gekeimt’ hat im Modell weniger Gewicht als ‘10 von 100’ Samen gekeimt. Stattdessen verwenden wir die Funktion ‘cbind()’ :
Dann können wir die Funktion glm() wie bekannt nutzen. Bei Anteilsdaten nutzen wir als family ‘binomial’
##
## Call:
## glm(formula = Keimungsrate ~ Substrat * Art, family = binomial,
## data = Ansaaten)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.61803 0.09285 17.427 < 2e-16 ***
## SubstratSteinwolle -1.57041 0.11570 -13.574 < 2e-16 ***
## ArtRadies 1.61079 0.20275 7.945 1.95e-15 ***
## SubstratSteinwolle:ArtRadies 1.33731 0.26856 4.980 6.37e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1266.76 on 39 degrees of freedom
## Residual deviance: 561.64 on 36 degrees of freedom
## AIC: 705.33
##
## Number of Fisher Scoring iterations: 5
Die summary des Modells zeigt, dass die Daten overdispersed sind. Auch hier gibt es für diesen Fall die Möglichkeit als family ‘quasibinomial’ anzugeben:
Zunächst lassen wir die Interaktion zwischen Substrat und Art weg, indem wir das Multiplikationszeichen durch ein Plus ersetzen
Jetzt vergleichen wir die Modelle mittels anaova(). Als Test muss jetzt der F-Test angegeben werden
## Analysis of Deviance Table
##
## Model 1: Keimungsrate ~ Substrat * Art
## Model 2: Keimungsrate ~ Substrat + Art
## Resid. Df Resid. Dev Df Deviance F Pr(>F)
## 1 36 561.64
## 2 37 584.97 -1 -23.326 1.603 0.2136
Wie wir sehen, ist die Interaktion nicht signifikant (p > 0.05) Also vereinfachen wir das Modell weiter und nehmen jetzt das Substrat als erklärende Variable raus
Modell3 <- glm(Keimungsrate ~ Art, family = quasibinomial, data = Ansaaten)
# Wieder vergleichen:
anova(Modell2, Modell3, test = "F")## Analysis of Deviance Table
##
## Model 1: Keimungsrate ~ Substrat + Art
## Model 2: Keimungsrate ~ Art
## Resid. Df Resid. Dev Df Deviance F Pr(>F)
## 1 37 584.97
## 2 38 767.96 -1 -182.99 11.644 0.001573 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Das Substrat hat einen hochsignifikanten Effekt auf die Keimungsrate. Als letztes testen wir noch die Art
Jetzt müssen wir Modell4 natürlich gegen Modell2 testen. Die Modelle müssen immer ‘genestet’ sein, sonst kann kein Signifikanzniveau berechnet werden. Genestet bedeutet, dass ein Modell eine einfachere Variante (eine oder mehrere Variable/n oder Interaktion/en weniger) des anderen Modells ist.
## Analysis of Deviance Table
##
## Model 1: Keimungsrate ~ Substrat + Art
## Model 2: Keimungsrate ~ Substrat
## Resid. Df Resid. Dev Df Deviance F Pr(>F)
## 1 37 584.97
## 2 38 1108.34 -1 -523.38 33.303 1.283e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Auch die Art hat einen hochsignifikanten Einfluss auf die Keimungsrate. Damit ist die Analyse der Keimungsraten abgeschlossen.
Bonitur-Daten sind recht speziell, weil es keine metrischen Daten sind, sondern Ordinal-Daten und darüber hinaus auch einigermaßen subjektiv. Trotzdem sind hierfür Methoden zur Auswertung entwickelt worden. Sie stammen aus der Psychologie, in der häufig Untersuchungen stattfinden, bei denen Personen zum Beispiel ihre Zufriedenheit auf einer Skala von 1 bis 5 bewerten sollen.
Wir schauen uns dazu exemplarisch die Bonitur-Daten zur Entwicklung der Bohnen an. Dazu filtern wir den Datensatz Ansaaten, indem wir nur die Zeilen nehmen, in denen das Argument Art == “Bohne” zutrifft (das doppelte Gleichzeichen wird immer benutzt, wenn eine Abfrage gemacht wird, ein einfaches Gleichzeichen ist in R eine Zuweisung also synonym mit dem Pfeil <- ). Mit dem Pipe-Operator %>% verbinden wir diese Auswahl mit der Auswahl der Spalten, die wir mit der Funktion select() machen. Hier können wieder einfach die header (=Spaltenüberschriften) der benötigten Spalten angegeben werden. Vorher laden wir noch zwei Pakete, die für die Auswertung von Ordinal-Daten entwickelt wurden
## Loading required package: xtable
##
## Attaching package: 'likert'
## The following object is masked from 'package:car':
##
## recode
## The following object is masked from 'package:dplyr':
##
## recode
##
## Attaching package: 'ordinal'
## The following object is masked from 'package:dplyr':
##
## slice
Bohne_Entw_l <- filter(Ansaaten, Art == "Bohne") %>% select(Substrat, ID, Entwicklung)
head(Bohne_Entw_l)## Substrat ID Entwicklung
## 1 Presstopf 1 8
## 2 Presstopf 2 5
## 3 Presstopf 3 8
## 4 Presstopf 4 8
## 5 Presstopf 5 4
## 6 Presstopf 6 5
Als erstes sollten die Daten wieder geplottet werden. Bei Bonitur-Daten eignen sich Histogramme am besten. Um das zu tun, müssen die Bonitur-Noten in Faktoren umgewandelt werden. Für die weitere Analyse geben wir auch direkt an, wieviele mögliche Levels (verschiedene Noten) es gibt und dass sie eine aussagekräftige Reihenfolge haben (ordered = TRUE):
Jetzt können die Histogramme erstellt werden:
ggplot(Bohne_Entw_l, aes(Entwicklung)) +
geom_bar() +
facet_wrap(~Substrat, ncol=1) +
xlab("Bonitur Note") +
ylab("Anzahl")Das package ‘likert’ stellt einige Plot-Funktionen zur Verfügung, die speziell für Bonitur- und andere Typen von likert-Daten geeignet sind. Dazu brauchen wir die Daten allerdings im ‘wide’ Format. Im Moment haben wir die Daten im ‘long’ Format: es gibt eine eigene Spalte für das Substrat, sodass ‘Presstopf’ und ‘Steinwolle’ je 10 Mal genannt werden. Wir nutzen wieder die Funktion ‘reshape’ wie oben, jetzt allerdings in umgekehrter Richtung, vom ‘long’ zum ‘wide’ Format:
Bohne_Entw_w <- reshape(Bohne_Entw_l, v.names="Entwicklung", idvar = "ID",
timevar = "Substrat", direction="wide")
head(Bohne_Entw_w)## ID Entwicklung.Presstopf Entwicklung.Steinwolle
## 1 1 8 5
## 2 2 5 5
## 3 3 8 5
## 4 4 8 3
## 5 5 4 2
## 6 6 5 1
‘idvar’ sind jetzt also die Zeilenbezeichnungen, ‘timevar’ die Spaltenbezeichnungen, ‘v.names’ die eigentlichen Werte.
Als nächstes müssen wir ein sogenanntes ‘likert-Objekt’ aus den beiden Spalten mit den Bonitur-Noten erstellen, dann kann die plot-Funktion des likert-packages genutzt werden:
Dies ist eine andere - etwas kompaktere - Darstellung der gleichen Daten. Die %-Werte geben an, wieviel % im schlechten Bereich (1-3), im mittleren Bereich (4-6) und im guten Bereich (7 - 9) waren.
Bonitur-Daten sind Ordinal-Daten und haben damit andere Eigenschaften als metrische Daten. Wie oben erwähnt, erfüllen sie häufig nicht die Annahmen, die für eine normale Varianzanalyse zutreffen müssen. Wir greifen hier auf ‘cumulative link models’ (CLM) aus der library ‘ordered’ zurück. CLMs sind im Prinzip das gleiche wie proportional odds model, falls Ihnen dieser Begriff mal über den Weg läuft. Ähnlich wie in den vorherigen beiden Kapiteln formulieren wir wieder ein komplexeres Modell mit Substrat als erklärender Variable und das Null-Modell, indem angenommen wird, dass alle Bohnen sich ungeachtet der Behandlung gleich Entwickeln:
Model1 <- clm(Entwicklung ~ Substrat, data = Bohne_Entw_l, threshold = "equidistant")
Model2 <- clm(Entwicklung ~ 1, data = Bohne_Entw_l, threshold = "equidistant")Als threshhold geben wir ‘equidistant’ an, weil wir davon ausgehen, dass zum Beispiel der Unterschied in der Entwicklung zwischen den Bonitur-Noten 2 und 3 genauso groß ist, wie der Unterschied zwischen den Noten 3 und 4 usw.
Auch für CLMs gibt es die Funktion anova(). Hier wird allerdings kein t- oder F-Test durchgeführt, sondern ein Chi-square Test, der uns einen p-Wert liefert:
## Likelihood ratio tests of cumulative link models:
##
## formula: link: threshold:
## Model2 Entwicklung ~ 1 logit equidistant
## Model1 Entwicklung ~ Substrat logit equidistant
##
## no.par AIC logLik LR.stat df Pr(>Chisq)
## Model2 2 78.001 -37.000
## Model1 3 73.330 -33.665 6.6711 1 0.009799 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Wie sie sehen, ist der p-Wert etwa 0.01, das heißt, das Substrat hat einen signifikanten Effekt auf die Entwicklung der Bohnen. Hiermit ist auch eine grundlegende Analyse der Bonitur-Daten abgeschlossen.
Ein interessanter Link zum Thema, inklusive Auswertungsmöglichkeiten mit R: ‘Do not use averages with Likert scale data’
Im Gartenbau und in den Agrarwissenschaften allgemein haben wir es
viel mit Wachstum zu tun: Pflanzen wachsen, Früchte wachsen, Nutztiere
wachsen. Auf dem Campus Klein Altendorf beschäftigen sich
Wissenschaftler zum Beispiel mit dem Wachstum von Äpfeln und seine
Abhängigkeit von verschiedenen Umweltfaktoren. Dazu gibt es unter
folgendem Link ein Video.
In diesem Kapitel schauen wir uns an, wie solche Zeitreihen mit R
dargestellt werden können und wie man Unterschiede im Wachstum (oder
anderen über die Zeit aufgenommenen Parametern) analysieren kann. Am
Ende des Kapitels
Zunächst laden wir die library ‘ggplot2’, lesen die Wachstumsdaten aus der Datei ‘Braeburn.csv’ ein (im Order ‘data’ auf eCampus) und kontrollieren den Import. ID und Zeitpunkt wandeln wir in Faktoren um. Wie Sie sehen, haben wir Daten von 4 Früchten, 2 von Ästen mit niedrigen Behangdichten und 2 von Ästen mit hohen Behangdichten, deren Umfang stündlich gemessen wurde.
## Datum Umfang ID Behang
## Length:2508 Min. : 1434 Min. :1.00 Length:2508
## Class :character 1st Qu.: 5049 1st Qu.:1.75 Class :character
## Mode :character Median : 7294 Median :2.50 Mode :character
## Mean : 7443 Mean :2.50
## 3rd Qu.: 9583 3rd Qu.:3.25
## Max. :14609 Max. :4.00
## Zeitpunkt
## Min. : 1
## 1st Qu.:157
## Median :314
## Mean :314
## 3rd Qu.:471
## Max. :627
Um das Wachstum dieser 4 Früchte zu visualisieren, nutzen wir die ggplot Funktion und geben bei den Aestetics (ein Parameter der Funktion) den Zeitpunkt (x-Achse), den Umfang dividiert durch 1000 um von Mikro- auf Millimeter zu kommen (y-Achse) und den Behang (niedrig/hoch) für die Gruppierung durch Farbe an:
ggplot(Braeburn, aes(Zeitpunkt, Umfang/1000, color=Behang)) +
geom_point() + labs(y="Umfang (mm)", x="Zeit") +
theme_bw()Im besten Fall, insbesondere wenn wir Unterschiede im Wachstum durch verschiedene Behandlung statistisch nachweisen wollen, haben wir nicht nur 2 Wiederholungen wie hier sondern deutlich mehr. Dies ist im Datensatz ‘Hühner’ der Fall, in dem das Gewicht von Junghennen, die mit unterschiedlichem Futter aufgezogen wurden, über 21 Tage protokolliert wurde.
Auch diesen Datensatz lesen wir ein (auch Ordner ‘data) und wandeln ’Futter’ und ‘ID’ in Faktoren um:
## X Gewicht Tag ID Futter
## 1 1 42 0 1 1
## 2 2 51 2 1 1
## 3 3 59 4 1 1
## 4 4 64 6 1 1
## 5 5 76 8 1 1
## 6 6 93 10 1 1
## X Gewicht Tag ID
## Min. : 1.0 Min. : 35.0 Min. : 0.00 Min. : 1.00
## 1st Qu.:145.2 1st Qu.: 63.0 1st Qu.: 4.00 1st Qu.:13.00
## Median :289.5 Median :103.0 Median :10.00 Median :26.00
## Mean :289.5 Mean :121.8 Mean :10.72 Mean :25.75
## 3rd Qu.:433.8 3rd Qu.:163.8 3rd Qu.:16.00 3rd Qu.:38.00
## Max. :578.0 Max. :373.0 Max. :21.00 Max. :50.00
## Futter
## Min. :1.000
## 1st Qu.:1.000
## Median :2.000
## Mean :2.235
## 3rd Qu.:3.000
## Max. :4.000
Das entsprechende Diagramm sieht folgendermaßen aus:
ggplot(Hühner, aes(Tag, Gewicht, color=Futter)) +
geom_point() + labs(y="Gewicht (g)", x="Tage seit Geburt") +
theme_bw()Für eine bessere Übersicht wäre es wünschenswert, nicht alle Datenpunkte zu plotten¸ sondern den Mittelwert des Gewichts der Hühner einer Behandlung und ein Maß entweder für die Varianz in den Daten oder die Genauigkeit des Mittelwertes. Mit der ggplot-Funktion ‘stat_summary’ können wir genau das machen: der Wert ‘mean_se’ im Parameter ‘fun.data’ bewirkt, dass Mittelwert und Standardfehler des Mittelwertes berechnet und geplottet werden:
ggplot(Hühner, aes(Tag, Gewicht, color=Futter)) +
stat_summary(fun.data=mean_se, geom="pointrange") +
labs(y="Gewicht (g)", x="Zeit (Tage seit Geburt)")Wir sehen in diesem Plot einen deutlichen Trend, nämlich dass das Futter des Typs 3 zu schnellerer Gewichtszunahme führt, als die übrigen. Die nächsten beiden Abschnitte erklären, wie wir testen, ob dieser Trend tatsächlich signifikant ist.
Im Prinzip sollte man meinen, dass sich Gewichts- oder Wachstumsdaten gut eignen sollten, um eine Varianzanalyse durchführen zu können: sie sind kontinuierlich und wir erwarten im Prinzip eine Normalverteilung der Residuen. Allerdings ist eine grundlegende Annahme verletzt: die Unabhängigkeit der Daten. Wenn wir immer wieder den gleichen Apfel oder das gleiche Huhn messen, können wir davon ausgehen, dass die Werte einander ähnlicher sind, als wenn wir immer andere Indiviuen nehmen würden. Wir haben also keine echten unabhängigen Wiederholungen sondern Pseudoreplikationen. Andere Beispiele die zu Pseudoreplikationen führen sind Versuchsblöcke, Standorte, unterschiedliche Probennehmer, etc.
Wie können solche Faktoren adäquat im Modell berücksichtigt werden? Der klassische Ansatz war lange, sie als ‘normale’ Variablen, also als festen Effekt mit ins Modell aufzunehmen. Das verbraucht aber viele Freiheitsgrade, weil für jeden Level des Zufallseffekts (also jedes Huhn) ein Parameter geschätzt wird (z.B. das mittlere Gewicht jedes Huhns, siehe Abbildung unten). Hier gibt es eine recht gute Erklärung für Freiheitsgrade: https://welt-der-bwl.de/Freiheitsgrade. Je weniger Freiheitsgrade es gibt, desto schwieriger ist es in einer Varianzanalyse, tatsächlich existierende Unterschiede in den Behandlungsgruppen auch als solche zu erkennen (Fähigkeit, Unterschiede zu erkennen = Power der Analyse). Deshalb ist der bessere Ansatz, die Variable als Zufallseffekt einzubeziehen. Das beruht auf der Annahme, dass die Basis-Parameter für jedes Level dieser Variable (jedes Huhn) aus einer gemeinsamen Verteilung stammen. Geschätzt werden muss für das Modell nur die Breite, also Varianz dieser Verteilung. Auf diese Weise wird nur 1 Freiheitsgrad verwendet und die Power der Analyse bleibt höher.
Feste Effekte
Zufällige Effekte
Um zufällige Faktoren zu berücksichtigen, benötigen wir gemischte Modelle.
Ein Paket, das Funktionen für gemischte Modelle zur Verfügung stellt ist ‘lme4’. Dieses Paket laden wir jetzt (nachdem Sie es installiert haben).
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
Die Funktion lmer() fitted lineare gemischte Modelle. Die abhängige Variable und die festen erklärenden Variablen werden wie bekannt in Zusammenhang gebracht: In unserem Beispiel steht das Gewicht links von der Tilde und Tag und Futter werden mit einem Multiplikationszeichen verbunden, um auszudrücken, dass in dem Modell beide Faktoren und ihre Interaktion einen Einfluss auf das Gewicht haben können. Dahinter kommen in runden Klammern die Zufallseffekte: Hier gehen wir davon aus, dass es einen Zufallseffekt geben könnte (individuelle Unterschiede im Gewicht der Hühner), der sich auf den y-Achsenabschnitt der Zeitreihen auswirkt. Das wird ausgedrückt, indem die ID hinter den senkrechten Strich geschrieben wird. Um die Modelle im Anschluss vergleichen zu können, müssen wir noch die Standardmethode ‘REML’ auf ‘FALSE’ setzen.
Danach geht es wie gewohnt weiter: wir vereinfachen das Modell, zuerst durch Ersetzen des Multiplikationszeichens mit einem Additionszeichen, womit wir die Interaktion zwischen Tag und Futter rauswerfen. Ist die Interaktion nicht signifikant, nehmen wir als nächstes das Futter ganz raus. Der Zeitpunkt (=Tag) sollte bei solchen Zeitreihen-Analysen allerdings immer im Modell bleiben, es sei denn, wir gehen davon aus, dass es möglicherweise gar keine Veränderung der abhängigen Variable über die Zeit gibt.
model2 <- lmer(Gewicht ~ Tag + Futter + (1 | ID), data=Hühner, REML = F)
model3 <- lmer(Gewicht ~ Tag + (1 | ID), data=Hühner, REML = F)Dann kann man die Modelle mittels anova() vergleichen. Wie Sie sehen, können auch mehrere Modelle gleichzeitig verglichen werden.
## Data: Hühner
## Models:
## model3: Gewicht ~ Tag + (1 | ID)
## model2: Gewicht ~ Tag + Futter + (1 | ID)
## model1: Gewicht ~ Tag * Futter + (1 | ID)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## model3 4 5630.3 5647.8 -2811.2 5622.3
## model2 7 5619.2 5649.7 -2802.6 5605.2 17.143 3 0.0006603 ***
## model1 10 5508.0 5551.6 -2744.0 5488.0 117.184 3 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sie sehen, dass sowohl der Einfluss des Futters auf den Gewichtsverlauf als auch die Interaktion zwischen Futter und Tag hochsignifikant sind. Die Vorhersagen des besten Modells (also model1) können wir nun noch zu dem Plot hinzufügen:
Nachdem wir in den letzten Kapiteln Varianzanalysen mit verschiedenen Datentypen besprochen haben, können Sie testen, ob eine erklärende Variable “signifikant” ist, also ihr Einbezug in ein Modell die übrigbleibenden Fehler (=Residuen) signifikant reduziert. Was Sie noch nicht wissen ist, welche der Level (=Behandlungsgruppen) der erklärenden Variable signifikant unterschiedlich sind. Zum Beispiel haben Sie beim Ansaaten-Versuch den Effekt von 4 verschiedene Substraten auf die Keimungsrate von 2 Arten getestet. Mit einer ANOVA finden Sie heraus, ob das Substrat grundsätzlich einen signifikaten Effekt auf die Keimungsrate hat. Bei WELCHEM Substrat oder welchen Substraten aber die Keimungsrate signifikant höher ist oder sind als bei den anderen, wissen Sie noch nicht. Um das herauszufinden, benutzen wir post-hoc Tests, also Tests, die nach der ANOVA durchgeführt werden.
Am Ende dieses Kapitels wissen Sie
Natürlich können wir jetzt jedes Substrat nochmal einzeln mit einer ANOVA gegen jedes andere Substrat testen und schauen, ob das Ergebnis siginifkant ist. Allerdings gibt es dabei ein Problem: jeder Test hat per Definition eine Irrtumswahrscheinlichkeit von 5% (wir nennen eine erklärende Variable signifikant, wenn der p-Wert, also die Wahrscheinlichkeit solche oder noch extremere Daten zu finden, obwohl die Null-Hypothese wahr ist, kleiner als 0.05 = 5% ist). Je mehr Vergleiche wir machen, desto höher ist also auch die Wahrscheinlichkeit, dass wir in der gesamten Auswertung eine Signifikanz finden, obwohl gar kein wahrer Effekt vorhanden ist - das Problem der multiplen Tests.
Einige Statistiker sagen deshalb, dass post-hoc Tests grundsätzlich vermieden werden sollten und auch nicht nötig sind, wenn das experimentelle Design des zugrundeliegenden Versuchs nur gut genug angelegt wurde. Die agrarwissenschaftlichen Realität sieht aber oft anders aus: es wird häufig eine Vielzahl von Behandlungslevels gegeneinander getestet und es interessiert uns nicht nur, OB die Wahl des Substrats einen Effekt auf die Keimungsrate einer Art hat, sondern auch bei welchem davon die Keimungsrate signifikant höher ist. Wahr ist aber, dass man gut durchdenken sollte, wann post-hoc Tests wirklich gebraucht werden und wie genau man sie durchführt.
Die generelle Herangehensweise ist, dass man für die Anzahl an Vergleichen, die man macht, korrigiert, indem die einzelnen Signifikanz-Niveaus so angepasst werden, dass man INSGESAMT auf eine Irrtumswahrscheinlichkeit von 5% kommt. Gängige Ansätze sind die Bonferroni-Korrektur und der Tukey’s HSD Test. Hier nutzen wir eine angepasste Bonferroni Korrektur, die sich ‘Holm’ nennt. Wer sich für die methodischen Details interessiert findet hier eine recht gute Erklärung.
Bitte laden Sie die Daten zum Ansaaten-Versuch von diesem Jahr
herunter: https://ecampus.uni-bonn.de/goto_ecampus_file_2859385_download.html
und installieren Sie das R-Paket ‘multcomp’.
## Loading required package: mvtnorm
## Loading required package: survival
## Loading required package: TH.data
## Loading required package: MASS
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
##
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
##
## geyser
library(tidyverse)
# Einlesen der Daten und Umwandeln der erklärenden Variablen in Faktoren
ansaaten <- read.csv("data/Ansaaten_neu.csv", header = TRUE)
ansaaten$Art <- as.factor(ansaaten$Art)
ansaaten$Substrat <- as.factor(ansaaten$Substrat)
# Kombinieren der beiden Spalten mit Anteilsdaten
ansaaten$keimung <- cbind(ansaaten$Angelaufen, ansaaten$Ausfälle)
# Umformen des Datensatzes, um einen Plot zu machen
Keimung_w <- reshape(ansaaten,
direction = "long",
varying = list(names(ansaaten)[4:5]),
v.names = "Anzahl",
idvar = c("Art", "Substrat", "Wdh"),
timevar = "Beobachtung",
times = c("Ausfaelle", "Angelaufen"))
head(Keimung_w)## Art Substrat Wdh keimung.1 keimung.2
## Bohne.Jiffy.1.Ausfaelle Bohne Jiffy 1 63 21
## Bohne.Jiffy.2.Ausfaelle Bohne Jiffy 2 34 50
## Bohne.Presstopf.1.Ausfaelle Bohne Presstopf 1 76 8
## Bohne.Presstopf.2.Ausfaelle Bohne Presstopf 2 76 8
## Bohne.Steinwolle.1.Ausfaelle Bohne Steinwolle 1 64 20
## Bohne.Steinwolle.2.Ausfaelle Bohne Steinwolle 2 70 14
## Beobachtung Anzahl
## Bohne.Jiffy.1.Ausfaelle Ausfaelle 63
## Bohne.Jiffy.2.Ausfaelle Ausfaelle 34
## Bohne.Presstopf.1.Ausfaelle Ausfaelle 76
## Bohne.Presstopf.2.Ausfaelle Ausfaelle 76
## Bohne.Steinwolle.1.Ausfaelle Ausfaelle 64
## Bohne.Steinwolle.2.Ausfaelle Ausfaelle 70
# Plot
ggplot(Keimung_w, aes(x=factor(Wdh), y=Anzahl, fill=Beobachtung)) +
facet_grid(. ~ Art * Substrat) +
geom_bar(position="stack", stat="identity") +
scale_fill_manual("legend",
values = c("Angelaufen" = "darkgreen", "Ausfaelle" = "lightgrey")) +
geom_col(position = position_stack(reverse = TRUE))# Generalisiertes lineares Modell unserer Hypothese:
# Substrat und Art haben einen Effekt auf die Keimungsrate
model1 <- glm(keimung ~ Substrat * Art, family=quasibinomial, data = ansaaten)
# Genestetes Modell ohne den Effekt von 'Substrat',
# um die Signifikanz dieser Variable testen zu können
model2 <- glm(keimung ~ Art, family = quasibinomial, data = ansaaten)
# Varianzanalyse der beiden Modelle
anova(model1, model2, test = "F")## Analysis of Deviance Table
##
## Model 1: keimung ~ Substrat * Art
## Model 2: keimung ~ Art
## Resid. Df Resid. Dev Df Deviance F Pr(>F)
## 1 8 31.78
## 2 14 326.94 -6 -295.16 12.977 0.0009685 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Mit der Funktion interaction() generieren wir jetzt eine
neue Spalte ‘gruppe’, in der die erklärenden Variablen kombiniert
werden. Mit diesen Behandlungsgruppen fitten wir ein neues glm, um sie
dann in einem nächsten Schritt in einer post-hoc Analyse gegeneinander
testen zu können.
## Art Substrat Wdh Angelaufen Ausfälle keimung.1 keimung.2 gruppe
## 1 Bohne Jiffy 1 63 21 63 21 Bohne.Jiffy
## 2 Bohne Jiffy 2 34 50 34 50 Bohne.Jiffy
## 3 Bohne Presstopf 1 76 8 76 8 Bohne.Presstopf
## 4 Bohne Presstopf 2 76 8 76 8 Bohne.Presstopf
## 5 Bohne Steinwolle 1 64 20 64 20 Bohne.Steinwolle
## 6 Bohne Steinwolle 2 70 14 70 14 Bohne.Steinwolle
Bevor wir das tun, sollten wir überlegen, welche Gruppen wir
tatsächlich vergleichen wollen. Die Abbildung oben legt nahe, dass Jiffy
stripes zu einer höheren Keimungsrate führen als Steinwolle und Sand und
wir wollen vielleicht testen, ob diese Unterschiede signifikant sind. In
der Funktion glht , die den post-hoc Test für uns
durchführt können wir festlegen, welche Gruppen wir gegeneinander testen
wollen:
summary(glht(model_posthoc,
linfct = mcp(gruppe =
#Ist die Differenz zwischen diesen Gruppen ungleich 0?
c("(Bohne.Jiffy)-(Bohne.Sand)=0",
"(Bohne.Steinwolle)-(Bohne.Sand)=0"))),test = adjusted("holm"))##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: User-defined Contrasts
##
##
## Fit: glm(formula = keimung ~ gruppe, family = quasibinomial, data = ansaaten)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## (Bohne.Jiffy) - (Bohne.Sand) == 0 -1.4319 0.5202 -2.753 0.0118 *
## (Bohne.Steinwolle) - (Bohne.Sand) == 0 -0.3725 0.5639 -0.661 0.5089
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- holm method)
Jetzt sehen Sie, dass der Unterschied zwischen Jiffy stripes und Sand bei der Bohne signifikant ist, der Unterschied Steinwolle - Sand aber nicht.
Signifikante Unterschiede, die zwischen den Behandlungsgruppen gefunden wurden, werden in Abbildungen häufig durch verschiedene Buchstaben gekennzeichnet. Zum Beispiel so (hier liegen andere Daten zugrunde):
Allerdings bietet R keine guten Standardlösungen dafür an. Wenn Sie den TukeyHSD Test durchführen, bietet das Packet ‘multcompView’ eine Möglichkeit. R-Code um die Abbildung oben zu erstellen finden Sie hier.
Auch im agrarwissenschaftlichen Bereich findet man zunehmend Studien, die nicht auf dem Formulieren von Null-Hypothesen und deren Testen mittels p-Werten beruhen (die frequentistische Statistik), sondern sich die Vorteile der Bayes-Statistik zunutze machen. Dazu gehören vor allem die intuitiver zu interpretierende Ergebnisse und die Möglichkeit, schon vorhandenes Wissen formalisiert in die Analyse einbeziehen zu können. Deshalb möchten wir Ihnen die Bayes-Statistik nicht vorenthalten! Das folgende Kapitel stützt sich stark auf einen sehr gut geschriebenen Artikel “Kurze Einführung in Bayes-Statistik mit R für Ornithologen” (2016) von Fränzi Korner-Nievergelt & Ommo Hüppop in ‘Die Vogelwarte’. Am Ende des Kapitels kennen Sie
Keine Einführung in die Bayesische Statistik kommt ohne den Satz von Bayes aus: Thomas Bayes war ein englischer Geistlicher, Philosoph und Statistiker von dem 1763 - posthum - ein Manuskript mit dem Titel ‘An essay towards solving a problem in the doctrine of chances’ veröffentlicht wurde. In diesem beschreibt er, wie basierend auf vorhandenem Wissen und zusätzlichen Beobachtungen X die Wahrscheinlichkeit dafür berechnet werden kann, dass eine Hypothese H zutrifft:
\[P(H|X) = \frac {P(X|H)×P(H)}{P(X)}\]
In Worten ausgedrückt: Die Wahrscheinlichkeit (P = probability) dass
die Hypothese zutrifft, nachdem wir die Daten betrachtet haben (P(H|X),
a-posteriori Wissen), ist gleich der
Wahrscheinlichkeit, die Daten so zu beobachten, angenommen die Hypothese
trifft zu (P(X|H), Likelihood), mal der
Wahrscheinlichkeit, dass die Hypothese zutrifft, bevor wir die Daten
betrachtet haben (P(H), a-priori Wissen oder
Prior), geteilt durch die Wahrscheinlichkeit die Daten
so zu beobachten ohne irgendeine Annahme zur Hypothese (P(X),
Normalisierungskonstante).
Wichtig ist noch anzumerken, dass Hypothesen dabei als Parameterwerte
ausgedrückt werden. Ist unsere Hypothese also, dass unter Behandlung A
der Mittelwert der abhängigen Variable größer ist als unter Behandlung
B, berechnen wir die Wahrscheinlichkeit für den Parameter Mittelwert(A)
- Mittelwert(B) und bekommen als Ergebnis, mit welcher
Wahrscheinlichkeit er > 0 ist.
In der Bayes-Statistik wird diese Idee verwendet, um bereits vorhandenes Wissen - das A-priori-Wissen (P(HP)) - mit den Information, die in den neu erhobenen Daten X stecken zu kombinieren und so das ‘geupdatete- A-posteriori-Wissen zu generieren. Wir sehen außerdem, dass wir als Ergebnis die Wahrscheinlichkeit unserer Hypothese bekommen. Das ist wesentlich leichter bekömmlich und näher an unserem ’normalen’ Denken als die Interpretation eines p-Wertes: ‘Die Wahrscheinlichkeit solche oder extremere Daten zu finden, wenn die Null-Hypothese wahr ist’.
Exkurs: Wahrscheinlichkeitsverteilungen
Um Wahrscheinlichkeiten darzustellen, zum Beispiel das a-priori Wissen, werden Wahrscheinlichkeitsverteilungen verwendet. In der foldenden Abbildung sehen Sie eine solche Wahrscheinlichkeitsverteilung: Je höher die Dichte, desto wahrscheinlicher der Parameterwert, der auf der x-Achse abgetragen ist. Die Gesamtfläche unter der Kurve beträgt immer 1 (genau irgendein Parameterwert trifft immer zu).
Wie berechnen wir also diese Wahrscheinlichkeit? Leider ist das meistens nicht ganz so trivial, wie es aussieht. Zwar können wir meistens die Likelihood berechnen (das wird auch in der frequentistischen Statistik zur Bestimmung der Parameterwerte unserer Modelle gemacht) und unser a-priori Wissen durch eine Wahrscheinlichkeitsverteilung festlegen. Schwierig wird es aber, zumindest sobald wir keine ganz simplen Modelle mehr haben, bei dem Teil P(X), der Wahrscheinlichkeit der Daten. Hierzu müsste man die Wahrscheinlichkeit der Daten über alle überhaupt möglichen Parameterwerte integrieren und das ist oft nicht möglich. Zum Glück kann dieses Problem aber mit einer Simulationstechnik, die in den 80er und 90er Jahren entwickelt wurde, umgangen werden:
MCMC-Verfahren können Wahrscheinlichkeitsverteilung annähern (in diesem Fall P(H|X)), die man nicht analystisch berechnen kann. Die intuitivste und ziemlich geniale Erklärung, die ich dazu gefunden habe ist die von Michael Franke auf seiner github Seite https://michael-franke.github.io/intro-data-analysis/Ch-03-03-estimation-algorithms.html Ich habe sie hier in Teilen übersetzt:
Wie die Äpfel an die Bäume kommen
Jedes Jahr im Frühling schickt Mutter Natur die Kinder los, um Äpfel in den Apfelbäumen zu verteilen. Für jeden Baum soll die Anzahl der Äpfel proportional zur Anzahl der Blätter sein: Riese Ralf mit doppelt so vielen Blättern wie Dünner Daniel soll also am Ende auch doppelt so viele Äpfel haben. Wenn es also insgesamt \(n_a\) Äpfel gibt, und \(L(t)\) die Anzahl der Blätter von Baum \(t\) ist, soll jeder Baum \(A(t)\) Äpfel bekommen.
\[A(t) = \frac {L(t)}{\sum L(t')}n_a\]
Das Problem ist nur, dass Mutter Natur \(\sum L(t')\) (die Normalisierungskonstante) nicht ausrechnen kann, also nicht weiß, wie viele Blätter alle Bäume zusammengenommen haben.
Die Kinder (Markov-Ketten) können aber zählen. Sie können zwar nicht alle Bäume besuchen und sich die Zahlen auch nicht so lange merken, aber immerhin für je zwei Bäume. Was sie jetzt tun, ist folgendes: sie starten je an einem zufälligen Baum (Parameterwert), hängen schonmal einen Apfel rein, zählen dort die Blätter \(L(t_1)\) und suchen dann nach einem Baum in der Nähe, von dem sie auch die Blätter zählen \(L(t_2)\) (der Zähler unserer Verteilung). Ist die Zahl der Blätter im zweiten Baum höher, gehen sie dort auf jeden Fall hin und hängen einen Apfel hinr ein. Ist sie niedriger gehen sie nur mit der Wahrscheinlichkeit \(L(t_2)/L(t_1)\) dorthin und hängen einen Apfel auf. So laufen sie weiter zwischen den Bäumen hin und her und am Ende sind die Äpfel ungefähr richtig verteilt. Je mehr Äpfel es sind, desto besser das Ergebnis. Die Besuchshäufigkeit der Kinder hat sich also der gewünschten Verteilung (die Mutter Natur nicht ausrechnen konnte) angenähert!
MCMC Verfahren machen das gleiche: eine MCMC Kette zieht zufällig
Werte für die Modell-Parameter und berechnet damit Ergebnis1 =
Likelihood der Daten * Prior. Dann zieht sie zufällig Werte, die um die
ersten Werte herum liegen, und berechnet dafür wiederum Likelihood der
Daten * Prior = Ergebnis2. Wenn Ergebnis2 höher ist als Ergebnis1,
springt sie dorthin und zieht von dort aus neue Parameter-Werte. Wenn
das Ergebnis2 niedriger ist, springt sie nur mit der Wahrscheinlichkeit
Ergebnis2/Ergebnis1 dorthin. Jetzt werden wieder zufällig Werte gezogen
usw.
In der folgenden Abbildung sind die erfolgreichen Sprünge durch blaue
Pfeile symbolisiert, die abgelehnten Sprünge durch grüne Pfeile.
Wenn dieser Prozess lang genug fortgesetzt wird, nähert sich die Verteilung der blauen Punkte der a-posteriori Wahrscheinlichkeitsverteilung der Parameter an, die wir haben wollen: Wir haben a-priori Wissen und die Information, die in den neu gesammelten Daten steckt, erfolgreich zum a-posteriori Wissen kombiniert!
Wir müssen uns entscheiden, welche Hühnerrasse wir in größerem Maßstab auf unserer Streuobstwiese halten möchten. Wir haben selbst schon 3 Kraienköppe und 6 Niederrheiner. Die Niederreihner sind uns ein bisschen lieber, weil sie weniger aggressiv sind, aber die Kraienköppe scheinen eine etwas höhere Legeleistung zu haben. Ist es deshalb wirklich sinnvoll, eher Kraienköppe zu wählen?
Schritt 1: Definition der A-priori-Wahrscheinlichkeitsverteilung
Wir erkundigen uns über die Legeleistung der beiden Rassen und erfahren, dass Niederrheiner zwischen 190 und 210 Eier pro Jahr legen, Kraienköppe zwischen 200 und 220. Entsprechend formulieren wir unsere a-priori Wahrscheinlichkeitsverteilungen:
## Loading 'brms' package (version 2.23.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').
# mit der Funktion set_prior definieren wir 2 Normalverteilungen, eine mit Mittelwert 200
# und Standardabweichung 1,5, und eine mit Mittelwert 210 und Standardabweichung 1,5.
# Unter 'coef' legen wir fest, zu welchen Parameter im Modell, das wir später
# fromulieren, diese Priors gehören.
priors <- c(set_prior("normal(200, 5)", coef = "rasseN"), set_prior("normal(210, 5)", coef = "rasseK"))
#So können wir sie plotten, um ein Gefühl dafür zu bekommen
curve(dnorm(x, 200, 5), from = 170, to = 230, xlab="Legeleistung", ylab="Dichte")2. Erhebung von neuen Daten
Jetzt geben wir unsere eigenen Beobachtungen der Legeleistung der 3 Kraienköppe und 6 Niederrheinern aus dem vorigen Jahr an:
rasse <- c("K", "K", "K", "N", "N", "N", "N", "N", "N")
ll <- c(225, 222, 221, 219, 219, 216, 221, 218, 217)
daten <- data.frame(rasse=rasse, ll=ll)3. Kombination der a-priori-Verteilung mit den Daten, um die a-posteriori-Verteilung zu erhalten
Um die a-posteriori Verteilung zu schätzen nutzen wir die Funktion ‘brm’ des Packets ‘brms’. Wie wir es schon von anderen Auswertungen kennen, formulieren wir zuerst unser Modell, nämlich das die Legeleistung ll von der Rasse abhängt. Die -1 in der Modell-Formel bewirkt, dass das Modell die Mittelwerte für die Legeleistung schätzt und nicht nur den Mittelwert für die erste Rasse und den Unterschied dazu. Unter ‘data’ geben wir unsere selbst erhobenen Daten ein und unter ‘prior’ die oben definierten Priors. Das argument ‘silent’ bestimmt, wie viele Informationen auf der Konsole ausgegeben werden, wenn die MCMC Ketten loslaufen. Diese Simulationen laufen je nach Modell-Komplexizität einige Sekunden oder Minuten.
4. Interpretation
Der obige Code hat die MCMC Simulation durchgeführt und wir können uns jetzt die a-posteriori Verteilungen für die Legeleistungen der Hühnerrassen anschauen:
## Family: gaussian
## Links: mu = identity
## Formula: ll ~ rasse - 1
## Data: daten (Number of observations: 9)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## rasseK 221.74 1.50 218.21 224.29 1.00 1991 1238
## rasseN 217.60 1.16 214.98 219.42 1.00 2180 1756
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 2.35 0.92 1.30 4.63 1.00 1313 1654
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
Die Zusammenfassung zeigt uns zuerst nochmal unsere Eingaben und ein paar Informationen zu den Markov-Ketten. Am interessantesten für uns sind natürlich die Schätzer für rasseK und rasseN und ihre Kredibilitätsintervalle. Wie Sie sehen, überschneiden sich die Konfidenzintervalle: der untere Wert von rasseK (l-95%) ist mit etwa 218 kleiner als der obere Wert (u-95%) von etwa 219 von rasseN, die legeschwächer ist.) Sigma ist die geschätzte Standardabweichung der angenommenen Normalverteilung und interessiert uns hier weniger.
In der Abbildung sind die Ergebnisse noch leichter zu interpretieren. Links sehen Sie die a-posteriori Verteilungen für beide Rassen (beachten Sie, dass die x-Achsen der beiden oberen Verteilungen nicht ausgerichtet sind). Rechts daneben die Dynamik der Markov-Ketten. Aus den a-posteriori Verteilungen ließe sich jetzt auch die genaue Wahrscheinlichkeit berechnen, dass die Legeleistung wirklich unterschiedlich ist. Da wir aber schon gesehen haben, dass die Kredibilitätsintervalle sich überschneiden, bleiben wir bei unserer Vorliebe für die Niederrheiner, auch wenn es nicht ausgeschlossen ist, dass sie ein paar Eier weniger legen.
| Frequentistisch | Bayesisch |
|---|---|
| Wahrscheinlichkeit für Daten, in Bezug auf Hypothese nur ja/nein Antwort | Wahrscheinlichkeit für Hypothesen |
| Vertrauensintervall: in welchem Intervall erwarten wir 95% der Mittelwerte von weiteren Zufallsstichproben gleicher Stichprobengröße | Kredibilitätsintervall: in welchem Wertebereich liegt mit einer Wahrscheinlichkeit von 95 % der wahre Mittelwert der Population |
| Annahme, dass es einen ‘wahren’ Wert gibt, die Beobachtungen aber zu einer Verteilung führen | Annahme, dass die beobachteten Werte wahr sind und eine intrinsische Varianz haben |
Die sinnvollere Interpretierbarkeit der Ergebnisse von bayesischen Analysen wird am besten durch folgene Abbildung dargestellt:
Bei einem unkritischen frequentistischen Ansatz würde das Ergebnis
lauten: “Der Effekt bei Gruppe/Behandlung A und E ist signifikant, also
Düngen wir jetzt mit dem getesteten Mittel. Bei den anderen Gruppen gibt
es anscheinend keinen Einfluss.”
Natürlich gibt es auch in B einen wichtigen Effekt und auch wenn der
Effekt in E statistisch signifikant ist, ist die Effektgröße auf den
Ertrag so klein, dass es von wirtschaftlichem Standpunkt aus
wahrscheinlich keinen Sinn macht, zusätzlich zu düngen. Die Ergebnisse
der bayesischen Analyse (also die Wahrscheinlichkeiten für die
Hypothese, dass es einen Effekt gibt), spiegeln diese Interpretationen
viel direkter wider.
Auch in diesem Kapitel gibt es eine Quiz-Fragen.